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Abstract 

Numerical studies of the two-dimensional Hubbard model have shown that it exhibits the 
basic phenomena seen in the cuprate materials. At half-filling one finds an antiferromagnetic 
Mott-Hubbard groundstate. When it is doped, a pseudogap appears and at low temperature 
d-wave pairing and striped states are seen. In addition, there is a delicate balance between 
these various phases. Here we review evidence for this and then discuss what numerical 
studies tell us about the structure of the interaction which is responsible for pairing in this 
model. 

1. Introduction 

A variety of numerical methods have been used to study Hubbard and t-J models and there 
are a number of excellent reviews. [1-8] The approaches have ranged from Lanczos diagonal- 
ization [1,2,9-11] of small clusters to density-matrix-renormalization-group studies of n-leg 
ladders [8,12-14] and quantum Monte Carlo simulations of two-dimensional lattices [3,15-24]. 
In addition, recent cluster generalizations of dynamic mean-field theory [4, 6, 7, 25-33] are 
providing new insight into the low temperature properties of these models. A significant 
finding of these numerical studies is that these basic models can exhibit antiferromagnetism, 
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stripes, pseudogap behavior, and d^^-y^ pairing. In addition, the numerical studies have 
shown how dehcately balanced these models are between nearly degenerate phases. Dop- 
ing away from half-filling can tip the balance from antiferromagnetism to a striped state in 
which half-filled domain walls separate 7r-phase-shifted antiferromagnetic regions. Altering 
the next-near- neighbor hopping t' or the strength of U can favor d^^^y^ pairing correlations 
over stripes. This delicate balance is also seen in the different results obtained using differ- 
ent numerical techniques for the same model. For example, density matrix renormalization 
group (DMRG) calculations for doped 8-leg t-J ladders find evidence for a striped ground 
state. [12] However, variational and Green's function Monte Carlo calculations for the doped 
t-J lattice, pioneered by Sorella and co-workers, [23,24] find groundstates characterized by 
(ij;2_j^2 superconducting order with only weak signs of stripes. Similarly, DMRG calculations 
for doped 6-leg Hubbard ladders [14] find stripes when the ratio of U to the near-neighbor 
hopping t is greater than 3, while various cluster calculations [27, 30-33] find evidence that 
antiferromagnetism and dx2_y2 superconductivity compete in this same parameter regime. 
These techniques represent present day state-of-the-art numerical approaches. The fact that 
they can give different results may reflect the influence of different boundary conditions or 
different aspect ratios of the lattices that were studied. The n-leg open boundary conditions 
in the DMRG calculations can favor stripe formation. Alternately, the cluster lattice sizes 
and boundary conditions can frustrate stripe formation. It is also possible that these differ- 
ences reflect subtle numerical biases in the different numerical methods. Nevertheless, these 
results taken together show that both the striped and the d^^^yi superconducting phases are 
nearly degenerate low energy states of the doped system. Determinantal quantum Monte 
Carlo calculations [21] as well as various cluster calculations show that the underdoped Hub- 
bard model also exhibits pseudogap phenomena. [27-32] The remarkable similarity of this 
behavior to the range of phenomena observed in the cuprates provides strong evidence that 
the Hubbard and t-J models indeed contain a signiflcant amount of the essential physics of 
the problem. [34] 

In this chapter, we will focus on the one-band Hubbard model. Section 2 provides an 
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overview of the numerical methods which were used to obtain the results that will be dis- 
cussed. We have selected three methods, determinantal quantum Monte Carlo, the dynamic 
cluster approximation and the density-matrix-renormalization group. In principle, these 
methods provide unbiased approaches which can be extrapolated to the bulk limit or in the 
case of the DMRG, to "infinite length" ladders. This choice has left out many other important 
techniques such as the zero variance extrapolation of projector Monte Carlo, [23,24], varia- 
tional cluster approximations, [25,26,29-32] renormalization group flow techniques, [35-37] 
high temperature series expansions [38-40] and the list undoubtedly goes on. It was moti- 
vated by the need to write about methods with which I have direct experience. 

In Section 3 we review the numerical evidence showing that the Hubbard model can in- 
deed exhibit antiferromagnetic, rf^2_j^2 pairing and striped low-lying states as well as pseudo- 
gap phenomena. From this we conclude that the Hubbard model provides a basic description 
of the cuprates, so that the next question is what is the interaction that leads to pairing in 
this model? Prom a numerical approach, this question is different from determining whether 
the groundstate is antiferromagnetic, striped, or superconducting. Here, one would like to 
understand the structure of the underlying effective interaction that leads to pairing. Al- 
though on the surface this might seem like a more difficult question to address numerically, it 
is in fact easier than determining the nature of the long-range order of the low temperature 
phase. The phase determination problem involves an extrapolation to an infinite lattice at 
low, or in the case of antiferromagnetism, to zero temperature. However, the pairing inter- 
action is short-ranged and is formed at a temperature above the superconducting transition 
so that it can be studied on smaller clusters and at higher temperatures. 

As discussed in Section 4, the effective pairing interaction is given by the irreducible 
particle-particle vertex F^p. Using Monte Carlo results for the one- and two-particle Green's 
functions, one can determine F^p and examine its momentum and Matsubara frequency 
dependence. [41,42] One can also determine if it is primarily mediated by a particle- hole 
magnetic {S — 1) exchange, a charge density {S — 0) exchange, or by a more complex 
mechanism. Section 5 contains a summary and our conclusions. 
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2. Numerical Techniques 



In this chapter we will be reviewing numerical results which have been obtained for the 2D 
Hubbard model. It would, of course, be simplest if one could say that these results were 
obtained by "exact" diagonalization on a sequence of LxL lattices with a "finite-size scaling" 
analysis used to determine the bulk limit. While one might not know the exact details of 
how this was done, one understands what it means. Unfortunately the exponential growth 
of the Hilbert space with lattice size limits this approach and other less familiar and often 
less transparent methods are required. 

In this chapter, we will discuss results obtained using the determinantal quantum Monte 
Carlo algorithm, [16, 43] a dynamic cluster approximation, [6, 27] and the density matrix 
renormalization group. [44,45, 5] All of these methods have been described in detail in the 
literature. However, to provide a context for the numerical results discussed in Sections 3 
and 4, we proceed with a brief overview of these techniques. 

Determinantal Quantum Monte Carlo 

The determinantal quantum Monte Carlo method was introduced in order to numerically 
calculate finite temperature expectation values. 

(A) = ill^ . (1) 
z 

Here, H includes —jiN so that Z — Tr e~^^ is the grand partition function. For the Hubbard 
model, the Hamiltonian is separated into a one-body term 

^ = XI {^l<r^j<^ + ^]c(^i<^ -fi^riia (2) 

(ij)a ia 

and an interaction term 

y = t/X(n,,-i) • (3) 

Then, dividing the imaginary time interval (0, (3) into M segments of width At, we have 
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In the last term, a Trotter breakup has been used to separate the non-commuting operators 
K and V. This leads to errors of order tUAr'^ which can be systematically treated by 
reducing the size of the time slice At. Then, on each — iAr shce and for each lattice site 

i, a discrete Hubbard- Stratonovi eh field [46] Si{Te) = ±1 is introduced so that the interaction 
can be written as a one-body term 

AtU 

^ Si{Te)=±l 

with A set by cosh(ATA) = exp(Ar^). The interacting electron problem has now been 
replaced by the problem of many electrons coupled to a r-dependent Hubbard-Stratonovich 

Ising field Si{T() which is to be averaged over. This average will be done by Monte Carlo 
importance sampling. 

For an L X L lattice, it is useful to introduce a one-body x lattice Hamiltonian 
ha{S{Ti)) for spin a electrons interacting with the Hubbard-Stratonovich field on the r^- 
imaginary time slice 

-/^X^^icT ± A^S'i(T£)nj(,. (6) 

i i 

The plus sign is for spin-up (a = 1) and the minus sign is for spin-down (cr = —1). Then, 
tracing out the fermion degrees of freedom, one obtains 

Z = det M| (S) det (S) . (7) 

{S} 

The sum is over all configurations of the (r^) field and M„[S) is an x matrix which 
depends upon this field, 

NUS) = I + Bl,B'l,_,---Bl. (8) 

/ is the unit x matrix and BJ — e'^'^^"^^^'^'-^^ acts as an imaginary time propagator 
which evolves a state from (-^ — 1) At to I At. 
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The expectation value (A) becomes 

iA}^j:A(S) detM|(S) detMi(S) 

{S} 

with A{S) the estimator for the operation A which depends only upon the Hubbard-Stratonovich 
field. Typically, we are interested in Green's functions. For example, the estimator for the 
one-electron Green's function is 

G^Are, (S)) = ^ ^ BJB!_, ■ ■ ■ Bl (10) 

and 

{S} 

The calculations of various susceptibilities and multiparticle Green's functions are straight- 
forward since once the Hubbard-Stratonovich transformation is introduced, one has a Wick 
theorem for the fermion operators. The only thing that one needs to remember is that dis- 
connected diagrams must be retained because they can become connected by the subsequent 
average over the 5*4 (r^) field. 

In computing the product of the B matrices, one must be careful to control round-off 
errors as the number of products becomes large at low temperatures or at large U where At 
must be small. In addition, there can be problems inverting the ill-conditioned sum of the 
unit matrix I and the product of the B matrices needed in Eqs. (jHl) and pO|) . Fortunately, 
a matrix stabilization procedure [16] overcomes these difficulties. 

For the half-filled case {fi = 0), provided there is only a near-neighbor hopping, H is 
invariant under the particle-hole transformation Cii = (— l)*c||. Under this transformation, 
the last term in Eq. © for a = —1 becomes 

-A5^5,(r,)(l-n,J (12) 

i 

SO that 

det Mi{S) = Y[ e^^"^"("^Met M^{S) (13) 
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This means that det M-^i^S) det Mi{S) is positive definite. In this case, the sum over the 
Hubbard-Stratonovich configurations can be done by Monte Carlo importance samphng, 
with the probabihty of a particular configuration {S{t£)} given by 

^^^^ ^ detMT(g)detM;(^) _ ^^^^ 

Given M- independent configurations, selected according to the probability distribution Eq. lTH 
the expectation value of A is 

{s} 

When the system is doped away from half-filling, the product detM^(S') detM|(S') can 
become negative. This is the so-called "fermion sign problem". In this case, one must use 
the absolute value of the product of determinants to have a positive definite probability 
distribution for the Hubbard-Stratonovich configurations. 

|detMK^)detM^(g)| 
"^^ E |detMt(^)detM|(^)| ^ ^ 

{S} 

Then, in order to obtain the correct results for physical observables, one must include the 
sign of the product of determinants [47] 

s = Sgn{detM^{S)detMi{S)) (17) 

in the measurement 

as) 

The II subscript denotes that the average is over configurations generated with the probability 
distribution Py given by Eq. (fT^ . If the average sign becomes small, there will be large 
statistical fiuctuations in the Monte Carlo results. For example, if = 0.1, one would 
have to sample of order 10^ times as many independent configurations in order to obtain the 
same statistical error as when {s)\\ = 1. On general grounds, one expects that decreases 
exponentially as the temperature is lowered. 

The average sign also decreases as U increases and makes it (exponentially) difficult 
to obtain results at low temperatures for U of order the bandwidth and dopings of interest. 
Figure H illustrates just how serious this problem is and why other methods are required. 
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Figure 1: The average of the sign of the product of the fermion determinants, Eq. El that 
enters in the determinantal Monte Carlo calculations is shown versus temperature for an 
8x8 lattice with U = 8t and (n) = 0.87. One can understand why calculations for U = 8t 
with T < 0.3 become extremely difficult. (Scalapino [19]) 

The Dynamic Cluster Approximation 

In the determinantal quantum Monte Carlo approach, one could imagine carrying out cal- 
culations on a set of LxL lattices and then scaling to the bulk thermodynamic limit. The 
dynamic cluster approximation [6] takes a different approach in which the bulk lattice is 
replaced by an effective cluster problem embedded in an external bath designed to represent 
the remaining degrees of freedom. In contrast to numerical studies of finite-sized systems 
in which the exact state of an L x L lattice is determined and then regarded as an ap- 
proximation to the bulk thermodynamic result, the cluster theories give approximate results 
for the bulk thermodynamic limit. Then, as the number of cluster sites increases, the bulk 
thermodynamic result is approached. 
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Figure 2: In the dynamic cluster approximation the Brillouin zone is divided into Nc cells 
each represented by a cluster momentum K. Then the self-energy and 4-point vertices are 
calculated on the cluster using an action determined by the inverse of the coarse-grained 
cluster-excluded propagator Eq. |2I1 This figure illustrates this coarse graining of the 
Brillouin zone for = 4: (Maier et al. [6]). 



In the dynamic cluster approximation, the Brillouin zone is divided into A^^^ = L"^ cells 
of size (27r/L)^. As illustrated in Fig. |2l each cell is represented by a cluster momentum K 
placed at the center of the cell. Then the self-energy S(k, is approximated by a coarse 
grained self-energy 

E(K + k',M;„) ~S,(K,cu„) (19) 
for each k' within the K*^ cell. The coarse grained Green's function is given by 

^(K, ^ t E w„-(.K-.k'-/.)-S.(K,.) ^2°) 

where the lattice self-energy is replaced by the coarse grained self-energy. Given G and S^, 
one can set up a quantum Monte Carlo algorithm [6,48] to calculate the cluster Green's 
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function. Here, the bulk lattice properties are encoded by using the cluster-excluded inverse 
Green's function 

g-\K, UJn) = G-\K, UJn) + S,(K, UJn) (21) 

to set up the bilinear part of the cluster action. In Eq. ()2H). the cluster self-energy has been 
removed from Q to avoid double counting. 
Then, the interaction on the cluster 

-j^ %+qtCktCk'_q| CK'i (22) 

K,K',Q 

is linearized by introducing a discrete r-dependent Hubbard-Stratonovich field on each r- 
slice and for each K point. In this way, the cluster problem is transformed into a problem of 
non-interacting electrons coupled to r-dependent Hubbard-Stratonovich fields. Integrating 
out the bilinear fermion field and using importance sampling to sum over the Hubbard- 
Stratonovich fields one evaluates the cluster Green's function Gc(K,a;„). From this, one 
evaluates the cluster self-energy 

S,(K, ^„) = g~\K, uJn) - G-\K, Un) . (23) 

Then, using this new result for Sc(K,co'„) in Eq. 1201 these steps are iterated to convergence. 

Measurements of correlation functions and the 4-point vertex are made in the same 
manner as for the determinantal Monte Carlo. That is, after the Hubbard-Stratonovich 
field has been introduced, one has a Wick's theorem for decomposing products of various 
time-ordered operators. However, in this case there is an additional coarse-graining of the 
Green's function intermediate state legs [6,42]. Since one is using a determinantal Monte 
Carlo method, there is also a sign problem for the doped Hubbard model. However, the 
self-consistent bath and the coarse-graining of the momentum space significantly reduce this 
problem so that lower temperatures can be reached. [49] 

The Density Matrix Renormalization Group 

The density-matrix-renormalization-group method was introduced by White. [44,45] Here, 
as illustrated in Fig. El for a one-dimensional chain of length L, the system under study is 
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Figure 3: The configuration of biocks used for tlie density matrix renormalization group 
algoritlim (Wliite [45]). 



separated into four pieces. A block Bi containing £ — L/2 — 1 sites on the left, a reflected 
Bf (right interchanged with left) block containing i' = L/2 — 1 sites on the right, and two 
additional sites in the middle. The left hand block Bi and its near-neighbor site are taken 
to be the system to be studied, while the block B^ plus its adjacent site are considered 
to be the "environment". The entire system is diagonalized using a Lanczos or Davidson 
algorithm to obtain the ground state eigenvalue and eigenvector ■0o- Then, one constructs a 
reduced density matrix from ipo 

Pii' = Yl ^oij^oi'j ■ (24) 
j 

Here, ■0oij = (^lOlV'o) with \i) a basis state of the £ + 1 block and \j) a basis state of the 
£' + 1 "environment" block. Then the reduced density matrix is diagonalized and m 
eigenvectors, corresponding to the m largest eigenvalues are kept. The Hamiltonian for 
the left hand block and its added site B'g_^^ is now transformed to a reduced basis consisting 
of the m leading eigenstates of pn'. The right hand environment block iJ^+i is chosen to 
be a reflection of the system block including the added site. Finally, a superblock of size 
L + 2 is formed using ii/^+i, i^^+i and two new single sites. Open boundary conditions are 
used. Typically, several hundred eigenstates of the reduced density matrix are kept and 
thus, although the system has grown by two sites at each iteration, the number of total 
states remains fixed and one is able to continue to diagonalize the superblock. 

This infinite system method suffers because the groundstate wave function ipo continues 
to change as the lattice size increases. This can lead to convergence problems. Therefore, 
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in practice, a related algoritlim in which the length L is fixed has been developed. In this 
case, instead of trying to converge to the infinite system fixed point under iteration, one has 
a variational convergence to the ground state of a finite system. This finite chain algorithm 
is similar to the one we have discussed but in this case the total length L is kept fixed and 
the separation point between the system and the environment is moved back and forth until 
convergence is achieved. [45, 5] Following this, one can consider scaling L to infinity. 

In a sense, the density-matrix-renormalization-group method is a cluster theory. It em- 
beds a numerical renormalization procedure in a larger lattice in which an exact diagonaliza- 
tion is carried out. The division of the chain into the system and the environment is similar 
in spirit to the embedded cluster and Q^^. The use of the reduced density matrix, corre- 
sponding to the groundstate, to carry out the basis truncation provides an optimal focus on 
the low-lying states. 

An important aspect of this approach is how rapidly the eigenvalues of the reduced density 
matrix fall off. This determines how many m states one needs to obtain accurate results. 
Unfortunately, for the study of n-leg Hubbard ladders, this fall-off becomes significantly 
slower as n increases and many more states must be kept. In addition, it appears that the 
behavior of the pairfield-pairfield correlation function is particularly sensitive to the number 
of states m that are kept. A measure of the errors associated with the truncation in the 
number m of density matrix eigenstates that are kept, is given by the discarded weight 

D 

Wm= 5^ (25) 

i=m+l 

Here, D is the dimension of the density matrix and Wi is its i^^ eigenvalue. A useful approach 
is to increase m and extrapolate quantities to their values as Wm ~^ 0. The error in the 
groundstate eigenvalue varies as Wm and a typical extrapolation is shown in Fig. |3] For 
observables which do not commute with H, the errors vary as y/Wm- 
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Figure 4: DMRG results for the ground state energy of a 2000-site Heisenberg spin-one-half 
chain versus the discarded weight Wm- The exact Bethe ansatz energy is shown as the line 
at the bottom of the figure. (S.R. White) 

3. Properties of the 2D Hubbard Model 

As we have discussed, the particle-hole symmetry of the half-filled Hubbard model with a 
near-neighbor hopping t leads to an absence of the fermion sign problem. In this case, the 
determinantal Monte Carlo algorithm [43] provides a powerful numerical tool for studying 
the low temperature properties of this model. In a seminal paper, Hirsch [15] presented 
numerical evidence that the groundstate of the half-filled two-dimensional Hubbard model 
with a near-neighbor hopping t and a repulsive on-site interaction U > had long-range 
antiferromagnetic order. In this work, simulations on a set of LxL lattices were carried out. 
For each lattice, simulations were run at successively lower temperatures and extrapolated 
to the T = limit. Then, a finite-size scaling analysis was used to extrapolate to the bulk 
T = limit. 

This work set the standard for what one would like to do in numerical studies of the 
Hubbard model. Unfortunately, the fermion sign problem prevents one from carrying out 
a similar determinantal Monte Carlo analysis for the doped case. However, various other 
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methods have been developed which provide information on tlie doped Hubbard model. 
Here, we will discuss results obtained from a dynamic cluster Monte Carlo algorithm. [6] This 
method also provides a systematic approach to the bulk limit as the cluster size increases. As 
noted in Sec. 2, the dynamic cluster Monte Carlo still suffers from the fermion sign problem, 
although to much less of a degree than the standard determinantal Monte Carlo. Maier et 
al. [27] have made the important step of studying the doped system on a sequence of different- 
sized clusters ranging up to 26 sites in size. Furthermore, this work recognized the importance 
of cluster geometry and developed a Betts'-like [50] grading scheme for determining which 
cluster geometries are the most useful in determining the finite-size scaling extrapolation. 

Following this, we review a density-matrix-renormalization-group (DMRG) study [14] of 
a doped 6-leg Hubbard ladder in which the authors extrapolate their results to the limit of 
zero discarded weight and to legs of infinite length. This work provides evidence that static 
stripes exist in the ground state for large values of U{U ^ 12t) but are absent at weaker 
coupling {U < 3t). We conclude this section with a discussion of the pseudogap behavior 
which is observed in the lightly doped Hubbard model when U is of order the bandwidth. 

The Antiferromagnetic Phase 

Determinantal quantum Monte Carlo results for the equal-time magnetization-magnetization 
correlation function 



with mf = — njj) are plotted in Fig. |3| These results are for a half- filled Hubbard 
model on a 10x10 lattice at a temperature T = O.lt with U = 4t. At this temperature, the 
antiferromagnetic correlation length exceeds the lattice size and the cluster is essentially in 
its groundstate. Strong antiferromagnetic correlations are clearly visible in C{i). 
The magnetic structure factor, shown in Fig. [Ul 



has a peak at g = (tt, tt) reflecting the antiferromagnetic correlations. As shown in Fig. [71 
as the temperature is lowered 3(7:, tt) grows and then saturates when the antiferromagnetic 
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Figure 5: The equal-time magnetization-magnetization correlation function C{ix,^y) on a 
10 X 10 lattice with U = 4t, (n) = 1 and T = O.lt. The horizontal axis traces out the 
triangular path on the lattice shown in the inset. Strong antiferromagnetic correlations are 
seen (Hirsch [15], Moreo et al. [17]). 



correlations extend across the lattice. If there is long-range antiferromagnetic order in the 
groundstate, the saturated value of ^(Tr, tt) will scale with the number of lattice sites = 
L X L. Furthermore, based upon spin-wave fluctuations, [51] one expects that the leading 
correction will vary as so that 



>g(7r,7r) 
hm — — — 



{m^)^ _A_ 
3 N-2 



(28) 



Fig. IHl shows 5'(7r,7r)/A^ versus A^^^ for U = At and one sees that the groundstate has 
long-range antiferromagnetic order. In his original paper, Hirsch [15] concluded that the 
groundstate of the half-filled 2D Hubbard model with a near-neighbor hopping t would have 
long-range antiferromagnetic order for f/ > 0. 
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Figure 6: S{q^,qy) versus g^, for {n) = 1, 
U = 4t and T = .167t. The solid line is a 
fit to guide the eye (Hirsch [15], Moreo et 
al. [17]). 
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Figure 7: The antiferromagnetic structure 
factor ^(vr, vr) for (n) = 1 and ?7 = 4 as 
a function of the inverse temperature (3 for 
various lattice sizes. ^(Tr, vr) saturates when 
the coherence length exceeds the lattice size 
(Hirsch [15], White et al. [16]). 



dx2-y2 Pairing 

The structure of the pairing correlations in the doped 2D Hubbard model were initially 
studied using the determinantal Monte Carlo method. The d-wave pairfield susceptibility 

P,= r dr {A,{t) Aim (29) 
Jo 

with 

= 77^ J2(-^y4,cU (30) 

^^^^ e,5 

was calculated. Here 5 sums over the four near-neighbor sites of i and (—1)"^ gives the H 1 — 

sign alteration characteristic of d-wave pairing. The doped Hubbard model has a fermion 
sign problem, so that the Hubbard-Stratonovich fields must be generated according to the 
probability distribution -P||(>S') given by Eq. (fT^ . In this case, it is essential to include the 
sign factor s in the evaluation of observables. The (red) circles in Fig. |H1 show results [47] 
for Pd{T) obtained on a 4x4 lattice with {n) = 0.875 and U = 4t. If the sign s is not 
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Figure 8: The zero-temperature limit of S'(7r,7r)/A^ versus l/N^^"^. The results extrapolate 
to a finite value as ^ oo implying that there is long range antiferromagnetic order in the 
groundstate of the infinite lattice (Hirsch [15], White et al. [16]). 

included, one obtains the (blue) squares. The neglect of this sign in early work [52] left the 
false impression that the Hubbard model did not support dx2^y2 pairing. 

As seen, when the sign is included, the d-wave pairfield susceptibility increases as the 
temperature is lowered. However, over the temperature range accessible to the determinantal 
Monte Carlo, it remains smaller than the U = result Pdo, shown as the (blue) dashed line 
in Fig. ^1 In Ref. [16], it was argued that this behavior was due to the renormalization of 
the single particle spectral weight and that the significant feature to note was that Pd{T) 
was enhanced over 

^^(^) = ^ 5Z ^(P) ^(-P) - • (31) 

pn 

Here, G{p) is the dressed single particle Green's function determined from the Monte Carlo 
simulation and corresponds to the contribution of a pair of dressed but non-interacting 
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Figure 9: The d-wave pairfield susceptibility 
Pd{T) (red circles) for a 4 x 4 lattice with 
U = At and (n) = 0.875 versus temperature 
T measured in units of the hopping t. The 
(blue squares) show the erroneous result that 
is found if the fermion sign is ignored. (Loh 
et al. [47]) 



Figure 10: The d-wave pair- field suscepti- 
bility Pd{T) is shown as the open (red) cir- 
cles. The open (green) squares show results 
for the "noninteracting" pair-field suscepti- 
bility Pd(T) calculated using dressed single- 
particle Green's functions, Eq. EH while the 
dashed (blue) curve is the noninteracting 
susceptibility P^o calculated with the bare 
Green's functions. (White et al. [16]) 



holes. The fact that Pd{T) lays below Pd{T) implies that there is an attractive (i^.2_j^2-pairing 
interaction between the holes. Pd{T) is shown as the (green) curve labeled with open squares 
in Fig. CHI 

In order to determine what happens at lower temperatures, Maier et al. [27] have deter- 
mined Pd{T) using a dynamic cluster approximation. In a systematic study, they provided 
evidence that the doped Hubbard model contained a pairing phase. In this work, 

the authors adapted a cluster selection criteria originally introduced by Betts et al. [50] in 
a numerical study of the 2D Heisenberg model. For the Heisenberg model, Betts et al. [50] 
showed that an important selection criteria for a cluster was the completeness of the "al- 
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Figure 11: Cluster sizes and geometries used by Maier et al. [27] in their study of the d-wave 
pair-field susceptibility. The shaded squares represent independent d-wave plaquettes within 
the clusters. In small clusters, the number of neighboring d-wave plaquettes listed in 
Table 1 is smaller than 4, i.e., than that for an infinite lattice. (Maier et al. [27]) 



lowed neighbor shells" compared to an infinite lattice. They found that a finite-size scaling 
analysis was greatly improved when clusters with the most complete shells were selected. 
For a d-wave order parameter, Maier et al. noted that one needs to take into account the 
non-local 4-site plaquette structure of the order parameter in applying this criteria. As illus- 
trated in Fig. ^2 the 4-site cluster encloses just one d-wave plaquette. Denoting the number 
of independent near-neighbor plaquettes on a given cluster by Z^, the 4-site cluster has no 
near-neighbors so that Za=0. In this case the embedding action does not contain any pair 
field fluctuations and hence Tc is over estimated. Alternatively, the 8A cluster has space for 
one more 4-site plaquette {Z^ = 1) and the same neighboring plaquette is adjacent to its 
partner on all four sides. In this case the phase fluctuations are over estimated and Tc is 
suppressed. For the 16B cluster, one has Z^i = 2 while Z^i = 3 for the oblique 16A cluster. 
Thus, one expects that the pairing correlations for the 16B cluster will be suppressed relative 
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to those for the 16A cluster. The number of independent neighboring d-wave plaquettes 
for the clusters shown in Fig. ^2 are listed in Table 1. 



Cluster 




Tjt 


4 


(MF) 


0.056 


8A 


1 


-0.006 


18A 


1 


-0.022 


12A 


2 


0.016 


16B 


2 


0.015 


16A 


3 


0.025 


20A 


4 


0.022 


24A 


4 


0.020 


26A 


4 


0.023 



Table 1: Number of independent neighboring d-wave plaquettes and the value T^. obtained 
from a linear fit of the pair-field susceptibility in Fig. ^l(Maier et al. [27]). 

Results for the inverse of the pair field susceptibility versus T for U = At and (n) = 0.9 
are shown in Fig. ^| As expected, the 4-site cluster results over estimate Tc and the results 
for the 8A and ISA clusters do not give a positive value for Tc. However, successive Zd = A 
results on larger lattices fall nearly on the same curve. These results suggest that the 2D 
Hubbard model with U = At and (n) = 0.9 has a d^-i-y-i pairing phase. The dynamic cluster 
approximation leads to a mean field behavior close to Tc. [53] Values of Tc obtained using a 
mean field linear fit of the low temperature data for the various clusters are listed in Table 1. 

If Tc ~ .02t and we take t = 0.2 eV, this gives Tc ~ 50K. We believe that Tc will 
increase with U, reaching a maximum when U is of order the bandwidth. In addition, 
we expect that the transition temperature is sensitive to the one-electron tight binding 
parameters. An example which illustrates this is known from density matrix renormalization- 
group calculations for the 2-leg Hubbard ladder. [54] Figure El shows an average of the 
rung-rung pairfield correlations 

D = J](A(. + £)At(0) (32) 
for a 2 X 16 Hubbard ladder versus the ratio of the rung to leg hopping parameters t±/t. 
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Figure 12: The inverse of the d-wave pair-field susceptibihty is plotted versus T/t for various 
clusters. Here U = At and (n) = 0.9. (Maier et al. [27]) 



Here 

= cJi|c|2j — c|ij^c|2| (33) 

creates a pair on the i^^ rung of the ladder. The pairing, as measured by D exhibits a 
maximum at a value of t±/t when the minimum of the antibonding band at fc^: = and 
the maximum of the bonding band at fc^; = vr approach the fermi surface. For the half-filled 
noninteracting system, this would occur when t±/t = 2. The doping and the interaction U 
leads to a reduction of this ratio and to a fiattening of the dispersion which further enhances 
the single particle spectral weight near the fermi energy. If one considers the antibonding 
band to have ky = n and the bonding band to have ky = 0, then this behavior is similar to 
increasing the single-particle spectral weight near (0, vr) and (vr, 0) in the 2D Hubbard model. 
One also sees that the largest peak in D occurs when U is of order the bandwidth. 

Having argued that the bandstructure plays a key role in determining T^, it is important to 
note that this raises a puzzle. State of the art LDA calculations, as Andersen and coworkers 
have shown, [55] can be folded down to give material specific near-neighbor t, next-near- 
neighbor t', etc. hopping parameters. For the one-band Hubbard model one would then 
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Figure 13: D versus t±/t for various values of C//i at a filling (n) — 0.9375 (Noack et al. [54]). 

have for the one-electron energy 

£k — — 2i(cos kx + cos ky) — At' cos k^ cos ky — . . . (34) 

Prom an analysis of a large number of hole-doped cuprates, it was found that Tc is correlated 
with the range of the intra-layer hopping. [56] For the one-band Hubbard model that we have 
discussed, this analysis imphes that should increase as t'/t becomes more negative. The 
opposite trend is seen in both dynamic cluster [57] and density matrix renormalization-group 
calculations. [13, 58] However, a projected fermion calculation [59] finds that t' enters the 
effective interaction and can lead to an increase in Tc which is consistent with the conclusions 
of Ref . [56] . The resolution of this puzzle represents an important open problem. 
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Stripes 

In a DMRG study of 7rx6 Hubbard ladders with 4r holes, Hager et al. [14] found that 
the ground state was striped for strong coupling values of U{U = 12t). Using a systematic 
extrapolation they gave evidence that such stripes exist in infinitely long 6-leg ladders. These 
studies also found that for small values of U{U = 3t) there were no stripes in the ground 
state. This work extended earlier work [60] on a 7x6 system with 4 holes which found that 
a well-defined stripe formed for f//t ~ 8 to 12. The absence of stripes for weak coupling is 
consistent with the fact that weak coupling renormalization group studies of the Hubbard 
model find no evidence of a stripe instability. [36, 37] 

Using the DMRG technique, the ground state expectation values of the hole-density 

6 

hix)=J2i^-{nix,y))) (35) 

y=l 

and the staggered spin density 

6 

s{x) =Y,i-'^r^'{M^,y) -ni{x,y)) (36) 

y=l 

were evaluated for 7rx6 ladders with 4r holes. Periodic boundary conditions were used for 
the 6-site direction and open boundaries were used in the leg direction. Results for the hole 
h{x) and spin s{x) densities for a 21x6 ladder doped with 12 holes are shown in Fig.^^ One 
sees from the modulation of the hole density along the leg direction that three stripes have 
formed. These stripes, each associated with 4 holes, run around the 6-site cylinder. In earlier 
t-J studies, [12] a preferred stripe density of half- filling was found and we believe that the 2/3 
filling seen in Fig. is a consequence of the restriction to 6 legs. Just as in the t-J ladder 
calculations, the staggered spin density undergoes a 7r-phase shift where the hole density is 
maximal. The finite staggered spin density is an artifact of the DMRG procedure in which 
no spin symmetry was imposed. This, along with the open boundary conditions which break 
the translational symmetry allow the charge and spin density structures to appear in the 
h{x) and s{x) expectation values. 
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Figure 14: The hole h{x) (circle) and staggered spin s{x) (square) densities in the leg x- 
direction are plotted for a 21 x 6 ladder with 12 holes for U = 12t (solid symbols) and 
U = 3t (open symbols). (Hager et al. [14]) 



While stripe-like structures are seen in Fig. for both U/t = 12 and U/t = 3, the 
amplitude of the charge density modulations for U/t = 3 are both smaller and less regular 
than the U/t = 12 results. As discussed in Sec. 2, DMRG results for operators which 
are non-diagonal in the energy basis are expected to deviate from their exact values by the 
square root of the discarded weight y/W^ as the number of basis states is increased. Thus, to 
determine whether there are stripes in the ground state of an infinite ladder, Hager et al. [14] 
extrapolated their results for a set of 7rx6 ladders to Wm — ^ and then took R = 7r ^ oo. 
They did this for the wave-vector power spectrum of the charge density 



with 



R 



— ^sm{k^x) {h{x)) . 



(37) 



(3^ 



For a ladder with a periodic array of stripes separated by 7 sites, the maximum contribution 
to is associated with the wave vector 

K 2r + 1 2 



TT 



R+1 7 



(39) 
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Figure 15: The amplitude of the power spectrum component \H{kD\/\/R for the hole density 
modulation. Results for a fixed number (6000 < m < 8000) of density-matrix eigenstates 
(squares) and results extrapolated to the limit Wm — ^ (circles) are shown as a function of 
the inverse ladder length 1/R for U = 12t (solid symbols) and U = 3t (open symbols). The 
dashed lines are linear fits. (Hager et al. [14]) 



and 

i^max = I H{kl) I OC /lo (40) 

as R goes to infinity. In Fig. HSl the amplitude H^^ixikl) / VR is plotted for U/t = 12 and 
U/t = 3 versus the inverse of the ladder length R^^. The solid squares show the results 
when a fixed number (6000 < m < 8000) of density-matrix eigenstates are retained. The 
solid circles are the extrapolated Wm —>■ {m —>■ oo) results for U/t = 12. Similar results are 
shown using open symbols for f//t = 3. When the Wm — > results are then extrapolated 
to R oo, one sees clear evidence for stripes when U/t = 12 and an absence of stripes for 
U/t = 3. Note the importance of the Wm —>■ extrapolation in determining the absence of 
stripes for [//t = 3. 

The Pseudogap 

Besides the antiferromagnetic, d-wave pairing, and striped phases, the cuprates exhibit a 
normal state pseudogap below a characteristic temperature T* when they are underdoped. 



26 



This pseudogap manifests itself in a variety of ways. [61] There is a decrease in the Knight 
shift, reflecting a decrease in the magnetic susceptibility. [62] This was interpreted in terms 
of the opening of a pseudogap in the spin degrees of freedom. Observations of a similar 
suppression in the tunneling density of states, [63] the c-axis optical conductivity [64] and 
the specific heat [65] made it clear that there was a pseudogap in both the spin and charge 
degrees of freedom. ARPES studies show that in the hole-doped materials, a pseudogap 
opens near the (vr, 0) antinodal regions while in the electron-doped materials, at the lowest 
dopings, it opens along the nodal direction near (|, f ). [66] The pseudogap appears in the 
underdoped region of the phase diagram and weakens as optimal doping is approached. If 
the Hubbard model is to contain the essential physics of the cuprates, it should exhibit the 
pseudogap phenomenon. 

Before looking for evidence of pseudogap behavior in the doped Hubbard model, it is 
useful to first look at the structure of the single particle spectral weight for the half-filled 
Hubbard model. An important paper on this was that of Preuss et al. [20] Here, determi- 
nantal Monte Carlo calculations of the finite temperature single particle Green's function 
G{k,T) were carried out on an 8x8 periodic lattice. The spectral weight 

A{k,u) = -- I'mG{k,iuJn ^ uj + iS) (41) 
vr 

was then determined using a numerical maximum entropy continuation. Results for the 
half-filled case with U = 8t and T = O.lt are shown in Fig. [TBk . Here, A{k,uj) is plotted 
versus u for various k values in the Brillouin zone. Figure ITM) summarizes these results using 
a standard "band structure" uj versus k plot in which the dark areas signify a large spectral 
weight. This work and related studies [67] showed that when U was of order the bandwidth 
or larger, there were four bands consisting of two incoherent upper and lower Hubbard bands 
and two quasiparticle-like, narrow bands nearer u = 0. The inner bands were found to have 
a dispersion set by J = At'^/U while the outer, upper and lower Hubbard bands, appear as 
an essentially dispersionless incoherent background. 

The left hand part of Fig. E| shows the single particle density of states for the half-filled 
case. Here, when the temperature is small compared to the exchange energy J ~ At'^/U, one 
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Figure 16: Single-particle spectral weight A{k,uj) for an 8 x 8 Hubbard model at half-filling 
(n) = 1 with U = 8t and T = O.lt. (a) A{k,uj) versus u for different values of k and (b) uj 
versus k plotted as a "band-structure" where sizable structure in A{k,u) is represented by 
the strongly shaded regions and peaks by error bars. (Preuss et al. [20]) 



clearly sees the broad upper and lower Hubbard bands and the narrow inner bands. When 
the system is hole-doped, the chemical potential moves down into the narrow coherent band 
that lays below = for the half-filled case and at the same time the upper coherent band 
loses weight and disappears as shown on the right hand side of Fig. El This is also seen in 
Fig. 1181 which shows A{k,uj) for (n) = .95 from Ref. [20]. Here, one sees that the dispersing 
band below u = in the insulator and the band that the holes are doped into as the system 
becomes metallic are quite similar. At the same time, the narrow dispersing band that lays 
just above u; = in the insulating state has lost most of its spectral weight. 
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Figure 17: On the left, the single particle density of states N{uj) versus uj for U = 8t and 
{n) = 1. On the right, N{iu) for the hole doped {n) = 0.875, U = 8t case at T = 0.33t 
(Scalapino [68]). 

For the doped system, the fermion sign problem limited the temperature to T = t/3 for 
the determinantal data shown in Fig. but later similar determinantal Quantum Monte 
Carlo runs at T = 0.25t and a filling of 0.95 found evidence for the formation of a pseudogap 
near (vr, 0). [21] In this work, the spin susceptibility was shown to have a large spectral 
weight at well defined spin excitations for the doping and temperature range in which the 
pseudogap appeared. There was no pseudogap in the overdoped {n) < 0.8 regions where the 
spectral weight of the spin susceptibility became broad and featureless. 

Dynamic cluster Monte Carlo calculations [28] with U = 6t and (n) = 0.95 find that 
the magnetic spin susceptibility exhibits a clear decrease below a temperature T = O.lt, as 
shown in the inset of Fig. and simulations at T = .06t gave the results for A{k,uj) shown 
in Fig. ^1 Here, a pseudogap has opened for k = (7r,0), while the nodal region with k = 

|) remains gapless. In addition, a variety of other cluster calculations [7,29,30,32] have 
found pseudogap behavior in both hole- and electron-doped Hubbard models and studied 
its dependence on the next near-neighbor hopping t'. The t' dependence as well as the 
doping dependence are consistent with renormalization-group calculations which show the 
importance of umklapp scattering processes [37] and the short range antiferromagnetic spin 
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Figure 18: The single-particle spectral weight A{k,u) for the hole doped (n) = 0.95 system 
at a temperature T = 0.33t. This plot is similar to Fig. and shows what happens as holes 
are doped into the Mott-Hubbard insulator. (Preuss et al. [20]) 



correlations. 

In the next section, we turn to a discussion of the effective pairing interaction. Specifically, 
the structure of the two-particle irreducible vertex and its associated d-wave eigenfunction 
are analyzed. 



4. The Structure of the Effective Pairing Interaction 

As discussed in Section 3, determinantal quantum Monte Carlo studies of the doped two- 
dimensional Hubbard model find that dx2_y2 pairing correlations develop as the temperature 
is lowered and a dynamic cluster quantum Monte Carlo calculation on Betts' clusters finds 
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Figure 19: The single particle spectral weight A{k, u) versus u for the antinodal k = (vr, 0) 
(red curve) and nodal k = (7r/2,7r/2) (black curve) momenta of an underdoped (n) = 0.95, 
U = 6t Hubbard model at T = 0.1 It. The inset shows the temperature dependence of the 
magnetic susceptibility. (M. Jarrell) 

evidence for a finite temperature d-wave superconducting phase. Here we discuss how one 
can use numerical techniques to determine the structure of the interaction responsible for 
the pairing. The basic idea is to numerically calculate the four-point vertex F and the 
single particle propagator G (solid lines) shown in Fig. 1201 Then, using the particle-particle 
Bethe-Salpeter equation fFig. 120b). one can extract the two-particle irreducible vertex F^p 
which is the pairing interaction. As we will discuss, the four-point vertex F also contains 
information on the particle-hole magnetic {S = 1) and charge {S = 0) channels. Thus, it 
provides a natural framework for understanding the relationship of the pairing interaction 
to these other channels. 

Using quantum Monte Carlo simulations, one can calculate both the one- and two-fermion 
Green's functions 



G{x2,Xi) 



(rc^(x2)4(xi)) 



(42) 



and 



G2{x4,X3,X2,Xi) = - {T Ca^{x4) Cas{x3) CI^{X2) cl^{Xi 



» 



(43) 
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Figure 20: Bethe-Salpeter equations for (a) the particle-particle and (b) the particle- hole 
channels showing the relationship between the full vertex, the particle-particle irreducible 
vertex F^^, and the particle-hole irreducible vertex F^'', respectively, (c) Decomposition of 
the irreducible particle-particle vertex F''^ into a fully irreducible two-fermion vertex Airr 
plus contributions from the particle-hole channels (Maier et al. [42]). 



Here, c}^{xi) creates an electron with spin a at site xi and imaginary time r^. T is the usual 
r-ordering operator and we have suppressed the a indices. Fourier transforming on both the 
space and imaginary time variables, one obtains G{p) and 

T 

+ j;^ ^Pi+P2,P3+P4G{pi)G{p3)T{p^,p3; p2,pi)G{p2)G{pi) (44) 

with p = (p, iuJn). Then, using the Monte Carlo results for G and G2, one can determine the 
four-point vertex F from Eq. (jl^ . 

Given F and G, one can solve the Bethe-Salpeter equations shown in Fig. I^Uk and b to 
obtain the irreducible particle-particle and particle-hole vertices F^p and T^^. For example, 
in the zero center of mass and energy channel, the particle-particle Bethe-Salpeter equation 
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shown in Fig. I^Uk gives 



T{p'\p) = rPP(p'b) - ^rPP(p'|fc)GT(A;)Gi(-A;)r(A;b) 



T 



(45) 



with r(p'|p) = r{p',—p'] p,—p). Given F and G, one can then determine the irreducible 
particle-particle vertex F^p. This procedure is essentially the opposite of what one does in 
the traditional diagrammatic approach. There, one introduces an approximation for the 
irreducible vertex Fpp and solves Eq. for F. Here, we use Monte Carlo results for F 
and G and solve Eq. ()45|) for F^p. The Monte Carlo results for F satisfy crossing symmetry 
and Fpp(p'Ip) obtained from Eq. ()45|) is the effective particle-particle interaction. There is no 
approximation except for the fact that a finite lattice is used and one has the usual statistical 
Monte Carlo errors (and the small systematic finite At errors which can be eliminated by 
extrapolating At — > 0). 

The dominant pairing response, at low temperatures, is found to occur in the even fre- 
quency drc2_y2 channel. Since this channel is even in both the relative frequency and mo- 
mentum, it must be a spin singlet. Note that there are also spin singlet pairing channels 
which are odd in the relative frequency and momentum. However, the pairing instability in 
the doped Hubbard model comes from the even frequency and even momentum part of the 
irreducible particle-particle vertex. 



Determinantal quantum Monte Carlo results [41] for T^^{p'\p) obtained from an 8x8 lattice 
with U = 4t and (n) = 0.87 are shown in Fig. El Here, FPP(p'|p) is plotted for various 
temperatures as a function of q = p' — p with p = (tt, 0) and u;„ = uJn' = ttT. One sees 
that as the temperature is lowered, F^p peaks at large momentum transfers. The size of the 
effective pairing interaction Fpp also depends upon the energy transfer Um = uJn' — (^n, and 
falls off with ujm on a scale set by the characteristic spin-fluctuation energy. 

To obtain a more intuitive picture of the way in which the local repulsive Uni-^nij^ Hubbard 
interaction can lead to an effective attractive pairing interaction in the singlet channel, it is 



FPP(p'b) = i [FPP(p'b) + FPP(-p'b)] . 



(46) 
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Figure 21: The even irreducible particle- 
particle vertex TPP{q, Um = 0) for q = p' — p 
and p = (tt, 0) versus momentum transfer q 
along the (1, 1) direction. Here U = 4t and 
(n) = 0.875. As the temperature decreases 
below the temperature where spin-spin cor- 
relations develop, the strength of the interac- 
tion is enhanced at large momentum trans- 
fers. (Bulut et al. [41]) 



Figure 22: The real-space structure of 
rPP(R) at a temperature T = 0.25t for U = 
At and (n) = 0.87. When the singlet elec- 
tron pair is separated by one lattice spac- 
ing, R = X or y, the interaction is attractive, 
while it is strongly repulsive when R = 
and the pair occupy the same site. (Bulut et 
al. [41]) 



useful to construct the real space Fourier transform 

rf(R) = ^ E (P , ^7rT■ p, ^vrT) . (47) 

p,p' 

Values for rPP(R) are shown in Fig. |221 with the distance R between the two fermions 
measured from the central point. If two fermions occupy the same site, spin- up and spin- 
down, rPP(R = 0) ~ 9.6t. That is, the effective pairing interaction is even more repulsive 
than the bare U = At onsite Coulomb interaction. However, if two fermions in a singlet state 
are on near-neighbor sites, the effective interaction rPP(R = x or y) ~ —0.5t is attractive. 

In order to determine the structure of the pairing correlations which are produced by 
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Fg , we turn to the homogenous Bethe-Salpeter equation 

-^T.^l'(p\p')G^(p')Gii-p')Mp') = KUp)- (48) 
p' 

The temperature dependence of the leading eigenvalue in the particle-particle channel is 
plotted versus the temperature in Fig. [2S1 When this eigenvalue reaches 1, it signals an 
instability into a superconducting phase. Here, U = 4t with (n) = 0.85 and we are showing 
results obtained using the dynamic cluster approximation [42] for the 24-site fc-cluster dis- 
cussed in Section 3. The distribution of k points for the 24-site cluster is shown in the inset 
of Fig. EHl Similar results for T > 0.25t have been obtained using the determinantal Monte 
Carlo algorithm on an 8x8 lattice. [41] 

The eigenfunction corresponding to the leading particle-particle eigenvalue is a singlet 
and its K dependence, plotted in the inset of Fig. OJ shows that it has dx^-y^ symmetry. 
The frequency dependence of this eigenfunction at the antinodal point K = (vr, 0) is shown 
in the main part of Fig. EH Here, (j){{iT , 0) , ujn) has been normalized so that at ujn = vrT 
its value is 1. It is even in Un as it must be for a d-wave singlet to satisfy the Pauli 
principle. Also shown in this figure is the cUm-dependence of the Q = (vr, tt) spin susceptibility 
x{Q,^^m) normalized by (x(Q)O) + x{Q ^ '^'^T)) / 2 for comparison with 0((7r, 0), tUn). The 
boson Matsubara frequency dependence, Um = ^itl-kT , of the susceptibility is seen to interlace 
with the fermion, ujn = {2n + l)7rT, dependence of the eigenfunction. The momentum and 
frequency dependence of 0rf^2_^2 (i^', tu) reflects the structure of the pairing interaction F^p. 
The numerical results show that F^p is an increasing function of momentum transfer and is 
characterized by a similar energy scale to that which enters the spin susceptibility x{Qi ^m)- 

In a similar way, one can use F and G to solve for the irreducible particle-hole vertex 
shown in Fig. EUb . The homogenous Bethe-Salpeter equation for the channel with center- 
of-mass momentum Q, Matsubara frequency ujm = and 2;-component of spin = is 

E + k- k' + Q, k') G^{k' + q) G^{k') (^Q^{k') = K{Q) (J^Qaik) . (49) 
k' 

The leading eigenvalue in the particle-hole channel occurs for Q = [it, tt) for the 24-site k- 
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Figure 23: Leading eigenvalues of the Bethe- 
Salpeter equations in various channels for 
U/t = 4 and a site occupation (n) = 0.85. 
The Q = (vr, tt), Um = 0, 5* = 1 mag- 
netic eigenvalue is seen to peak at low tem- 
peratures. The leading eigenvalue in the 
even singlet Q = (0,0), Um = particle- 
particle channel has d^-i-y-i symmetry and in- 
creases toward 1 at low temperatures. The 
largest charge density eigenvalue occurs in 
the Q = (0,0), ojm = channel and satu- 
rates at a small value. The inset shows the 
distribution of /c-points for the 24-site clus- 
ter. (Maier et al. [42]) 
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Figure 24: The Matsubara frequency de- 
pendence of the eigenfunction 4>d^2_ 2 (K, ojn) 
of the leading particle-particle eigenvalue 
of Fig. 1221 for K = (vr, 0) normalized 
to 0(K, vrT) (red). Here, Un = (2n + 
l)7rT with T = 0.125t. For compari- 
son, the Matsubara frequency dependence 
of the normalized magnetic spin suscepti- 
bility 2x(Q,a;^)/[x(Q,0) + x(Q,2vrT)] for 
Q = (vr, vr) versus Um = 2'rmTT is also shown 
(green). In the inset, the momentum de- 
pendence of the eigenfunction 4>d^2_ 2 i.^^ 
normalized to 0^ 2_ 2 ((0, vr), vrT) shows its 
dx'i-y-i symmetry. Here, Un = t^T and the 
momentum values correspond to values of K 
which lay along the dashed line shown in the 
inset of Fig. ESI (Maier et al. [42]) 



cluster and carries spin 1. Earlier determinantal quantum Monte Carlo studies [17] on 8x8 
lattices show that for this doping the peak response is, in fact, slightly shifted from (7r,7r), 
but the 24-site fc-cluster used in the dynamic cluster calculation lacks the resolution to show 
this. As seen in Fig. |2S1 for this doping, the antiferromagnetic eigenvalue initially grows 
as the temperature is reduced, peaking at low temperatures. The largest eigenvalue in the 
5* = charge density channel occurs for Q = (0, 0) and Um = 0. Its temperature dependence 
is also plotted in Fig. 1221 

Returning to the question of the structure of the irreducible particle-particle vertex F^p, 
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Figure 25: (a) The irreducible particle-particle vertex T^p versus q = K — K' for various 
temperatures with ujn = ^n' = ttT. Here, K = (vr, 0) and K' moves along the momentum 
values of the 24-site cluster which lay on the dashed line shown in the inset of Fig. 1221 Note 
that the interaction increases with the momentum transfer as expected for a d-wave pairing 
interaction, (b) The q-dependence of the fully irreducible two-fermion vertex Am- (c) The 
q-dependence of the charge density {S = 0) channel for the same set of temperatures, 
(d) The q-dependence of the magnetic {S = 1) channel |$m- (Maier et al. [42]) 



we have seen that F^p peaks at large momentum transfers and has a frequency dependence 
reflected in ,^{K,uJn) which is similar to the spin susceptibility. However, we would 

like to understand one further aspect. Is the dominant contribution to the d^^^-y-i pairing 
interaction associated with an S* = 1 particle-hole channel? Alternatively, for example, one 
could have a charge density 5 = channel or a more complicated multiparticle-hole exchange 
process such as that suggested by the spin-bag picture. [69] 

In order to address this, we will make use of the representation of F^p shown diagram- 
matically in Fig. I^Ub . Here, F^p is decomposed into a fully irreducible vertex Am plus 
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contribution from particle-hole exchange channels. Because of the spin rotation invariance 
of the Hubbard model, one can separate the particle-hole channels into a charge density 
S = contribution and a spin S = 1 magnetic part. For the even frequency and even 
momentum (singlet pairing) part of the irreducible particle-particle vertex, Eq. one has 

rf(p» = Ai„(p» + 1 Mp',p) + 1 ^Up',p) ■ (50) 

The subscripts d and m denote the charge density (5* = 0) and magnetic {S = 1) particle-hole 
channels respectively, with 

1 



^d/m{P - P; P, -P') - ^d^/miP - -P') 



+ ^d/m{p' + p; -p, -p') - r^/m(p' + p] -p, -p') 



(51) 



Here, on the right hand side, the center of mass and relative wave vectors and frequencies in 
these channels are labeled by the first, second and third arguments, respectively. 

Results for the irreducible particle-particle interaction F^p obtained from the 24-site dy- 
namic cluster approximation are shown in Fig. 123 As we have seen, when the temperature 
is lowered, F^p increases as the momentum transfer q = p' — p increases. Using the results 
for FP'^, F, and G one can calculate the contributions $rf from the 5 = charge density and 
$m for the 5 = 1 magnetic channels. Subtracting these from Fpp gives Ai^r and results for 
each of these contributions are shown in Fig. 123 The dominant d^^-yi pairing contribution 
to FPP clearly comes from the 5 = 1 channel. 

At larger values of U, 4-site fc-cluster calculations [70] of the temperature dependence 
of the dr^2_y2 eigenvalue for (n) = 0.85 and U = At, 8t and 12t are shown in Fig. ESI Over 
the temperature range [71] shown in Fig. |2ni the dx2_y2 eigenvalue is largest for U = 8t. 
This is consistent with the 2-leg ladder result shown in Fig. and the expectation that the 
maximum transition temperature occurs for U of order the bandwidth. The dx2-y2 eigen- 
function (pd^2_y2 (-^? ^n) has the expected d-wave K dependence and its Matsubara frequency 
dependence for U = 4t and 8t are shown in Fig. EZI Here, as before, we also show the 
Um dependence of the spin susceptibility xiQ^^m)- As U increases, both 0^ ^_ ^{K^ojn) 
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Figure 26: The dx2_y2-waYe eigenvalue 
versus temperature T/t for (n) = 0.85 with 
U = U (red), U = 8t (blue) and U = 12t 
(green). These results were obtained for a 
4-site /c-cluster (Maier et al. [70]). 
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Figure 27: The Matsubara frequency depen- 
dence of the (ia,2_y2 eigenfunction (j)d{K,u!n) 
with K = (vr, 0) for U = At, 8t and 12t. The 
spin susceptibility xiQy ^m) for Q = (vr, vr) is 
also shown for comparison. Here, 0^ and x 
are normalized in the same way as the results 
shown in Fig. 123 (Maier et al. [70]). 



and xiQj^m) fall off more rapidly, reflecting the reduction in the frequency scale set by 
J ^ AtyU. 

5. Conclusions 

The numerical studies of the Hubbard model that we have reviewed show that it exhibits 
the basic properties that are observed in the cuprate materials: antiferromagnetism, d^^^yi- 
pairing, stripes and pseudogap phenomena. Numerical methods have also been used to study 
the structure of the interaction responsible for pairing in the Hubbard model. As discussed in 
Section 4, this can be done by directly calculating the irreducible particle-particle vertex F^^ 
or by studying the momentum and frequency dependence of the gap function 0rf^2_^2 (i^', tu). 
The decomposition of F^^ showed that the dominant pairing interaction arose from a spin-one 
particle-hole exchange. The strength of F^^ was found to increase with momentum transfer 
leading to (ia,2„j^2-pairing. Alternately, the {cosKx — cos Ky) momentum dependence of the 
gap function 0rf^2_^2 (i^T, tu) and the similarity of its Un dependence to that of the Q = (vr, vr) 
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spin susceptibility leads to the same conclusion: the pairing interaction in the doped Hubbard 
model is repulsive on site, attractive between near-neighbor sites and retarded on a time scale 
set by the inverse of the spin-fluctuation spectrum. It is important to recognize that this 
spectrum includes a particle-hole continuum. 

Now, one can ask whether this interaction is actually the mechanism responsible for 
pairing in the high Tc cuprate materials and how one would know this from experiments? 
As far as the momentum dependence of the interaction is concerned, ARPES studies [66] 
of the /c-dependence of the gap along with a variety of transport [72] and phase dependent 
studies [73, 74] provide strong evidence for the nodal d-wave character of the gap. While it is 
known that the chains in YBCO lead to an admixture of s-wave [75-77] and the momentum 
regions probed are primarily along the fermi surface, there is good reason to believe from 
the observed /c-dependence of the gap that the pairing interaction is indeed repulsive on 
site and attractive for singlets formed between near neighbor sites. It will be interesting to 
compare calculations for an orthohrombic Hubbard model with experiments [76,77], to see if 
the observed /c-dependence of the gap can provide additional help in identifying the pairing 
mechanism. 

Another characteristic of the interaction is its frequency dependence. Here, less is known 
but it seems likely that the frequency dependence of the gap and renormalization parameter 
will provide important insight into the mechanism. As one knows, it was the frequency 
dependence of the gap for the traditional low superconductors that provided the ultimate 
flngerprint identifying the phonon exchange pairing interaction, although at the time few 
doubted that this was the mechanism. In the high Tc case, the initial hope was that the 
d-wave momentum dependence of the gap would provide a sufficiently precise fingerprint. 
However, this has not been the case. For example, the exchange of Big phonons is known 
to favor d-wave pairing, [55,78], although its overall contribution to is small within the 
standard theory. A two-band Cu-0 model, in which fluctuations in circulating currents 
provide a d-wave pairing mechanism, has also been proposed. [79] Even within the framework 
of the Hubbard model there are different views regarding the dynamics. In the "Plain- 
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Vanilla-RVB" picture, [80] it has been suggested that the dynamics is set by an energy scale 
associate with the Mott-Hubbard gap. [81] However, our numerical results support a picture 
in which the dominant contributions come from particle-hole excitations within the relatively 

narrow band that the doped holes enter giving an energy scale of several times J . While the 
spectrum of these excitations extends down to zero energy, the main strength is associated 
with a broad spin-fluctuation continuum. [82,83] Thus it seems likely that the dynamics will 
again be important in identifying the mechanism. 

In addition to the traditional electron tunneling [84] and infrared conductivity [85] mea- 
surements, ARPES experiments provide an important tool for probing the frequency depen- 
dence of the renormalization parameter and the gap. Advances in the energy and momentum 
resolution of both ARPES [66] and neutron scattering [86] along with material preparation 
techniques that allow ARPES and neutron scattering to be done on the same material are 
opening new opportunities. Various RPA-BCS approximations have been used to model both 
the ARPES [87, 88] and neutron scattering data. [89, 90] One would clearly like to extend 
the numerical Hubbard model studies so that they can be used in making such experimental 
comparisons. 

Finally, in addition to the frequency and momentum dependence of the interaction, there 
is the question of its strength. The estimate for the transition temperature in Sec. 3 with 
U = At was relatively small. As discussed, we believe that for larger values of U (of order 
the bandwidth) and a more optimal bandstructurc, will increase. Beyond this, the actual 
Cu-0 structure has additional exchange paths and it is known that t — J—U Hubbard ladders 
can exhibit stronger pairing correlations. [91] Nevertheless, the question of the strength of 
the pairing interaction remains. It is not that several times J isn't a wide spectral range 
compared to the phonon scale of the traditional low temperature superconductors or that the 
system isn't strongly coupled with U of order the bandwidth. Rather it is that the strong 
coupling has created a delicately balanced system. [92] As discussed in Sec. 3, different 
numerical methods on different lattices find evidence in one case for d-wave pairing and in 
another for stripes. Thus small changes in local parameters may alter the nature of the 
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correlations and there is a question regarding the role of inhomogeneity in the cuprates. 
An interesting theory of "dynamic inhomogeneity-induced pairing" is discussed in another 
chapter of this treatise. [93] In this approach, pairing from repulsive interactions appears 
as a mesoscopic effect and the phenomena of high temperature superconductivity is viewed 
as arising from the existence of mesoscale structures. [93, 94] Recent STM measurements of 
impurities and inhomogeneities in BSCCO are providing important new information on the 
question of the local modulation of the pairing and its strength. [95-97] 

Thus, two decades after Bednorz' and Miiller's [98] discovery of the high Tc cuprates the 
question of the pairing mechanism remains open. However, it is clear that the desire to 
understand these materials has driven dramatic advances in the experimental energy and 
momentum resolution of ARPES and neutron scattering and the energy and spatial resolu- 
tion of STM. It was also largely responsible for the development of a variety of numerical 
techniques which are providing new insights into the electronic properties of a wide class of 
strongly correlated materials. 

Acknowledgement 

I would like to acknowledge past graduate students N. Bulut, J.M. Byers, M. Jarrell, E. Loh, 
R. Melko, R.M. Noack, S. Quinlan, and R. Scalettar and postdocs F. Assaad, N.E. Bickers, 
L. Capriotti, E. Dagotto, T. Dahm, J. Preericks, M.E. Platte, R. Eye, J. Hirsch, M. Imada, 
P. Monthoux, A. Moreo, M. Salkola, H.B. Schuttler and S.R. White who have been wiUing 
to teach me new things for so long. I would also like to acknowledge the pleasure I have 
had discussing and working on this problem with A.V. Balatsky, W. Hanke, P. Hirschfeld, 
S.A. Kivelson, T. Maier and D. Poilblanc. Pinally I want to thank my long time collaborator 
R.L. Sugar for his insights and encouragement. This work was supported by NSP Grant 
DMR02-11166 and the Department of Energy under EG02-03ER46048. 



42 



Bibliography 



[1] 

[2] 
[3] 
[4] 
[5] 
[6] 
[7 



E. Dagotto, Rev. Mod. Phys. 66, 763 (1994). 

J. Jaklic and P. Prelovsek, Adv. Phys. 49, 1 (1999). 

N. Bulut, Adv. Phys. 51, 1587, (2002). 

A. Georges, G. Kotliar, W. Krauth and M.J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996). 
R.M. Noack and S.R. Manmana, AIP Conf. Proc. 789, 93 (2005), |cond-mat/051032l| 
T. Maier, M. Jarrell, T. Pruschke and M. Hettler, Rev. Mod. Phys. 77 1027 (2005). 
A. -M.S. Tremblay, B. Kyung and D. Senechal, cond- mat/051 1334[ 



E. Dagotto and T.M. Rice, Science 271, 618 (1996). 
[9] D. Poilblanc, Phys. Rev. B 48, 3368 (1993); Phys. Rev. B 49, 1477 (1994). 
[10] F. Becca, A. Parola and S. Sorella, Phys. Rev. B 61, R16287 (2000). 
[11] P.W. Leung, Phys. Rev. B 62, R6112 (2000). 
[12] S.R. White and D.J. Scalapino, Phys. Rev. Lett. 81, 3227 (1998). 
[13] S.R. White and D.J. Scalapino, Phys. Rev. B 60, R753 (1999). 

[14] G. Hager, G. Wellein, E. Jeckelmann and H. Fehske, Phys. Rev. B 71, 75108 (2005). 

[15] J.E. Hirsch, Phys. Rev. B 31, 4403 (1985). 

43 



44 



BIBLIOGRAPHY 



[16] S.R. White, D.J. Scalapino, R.L. Sugar, E.Y. Loh, J.E. Gubernatis and R.T. Scalettar, 
Phys. Rev. B 40, 506 (1989); S.R. White, D.J. Scalapino, R.L. Sugar, N.E. Bickers and 
R.T. Scalettar, Phys. Rev. B 39, 839 (1989). 

[17] A. Moreo, D.J. Scalapino, R.L. Sugar, S.R. White and N.E. Bickers, Phys. Rev. B 41, 
2313 (1990). 

[18] A. Moreo, Phys. Rev. B 48, 3380 (1993). 

[19] D.J. Scalapino, Modern Perspectives in Many-Body Physics, pp. 199-225, (Sixth Physics 
Summer School, Australian National University), ed. by M.P. Das and J. Mahanty, 
World Scientific (1994). 

[20] R. Preuss, W. Hanke and W. von der Linden, Phys. Rev. Lett. 75, 1344 (1995). 

[21] R. Preuss, W. Hanke, C. Grober and H.G. Evertz, Phys. Rev. Lett. 79, 1122 (1997). 

[22] D. Duff'y, A. Nazarenko, S. Haas, A. Moreo, J. Riera and E. Dagotto, Phys. Rev. B 56, 
5597 (1997). 

[23] S. Sorella, G. Martins, R. Bccca, C. Gazza, L. Capriotti, A. Parola and E. Dagotto, 
Phys. Rev. Lett. 88, 117002 (2002). 

[24] V.I. Anisimov, M.A. Korotin, LA. Nekrasov, Z.V. Pehelkina and S. Sorella, Phys. Rev. 
5 66, 100502 (2002). 

[25] M. Potthoff, Eur. Phys. J., B32, 429436 (2003). 

[26] C. Dahnken, M. Aichhorn, W. Hanke, E. Arrigoni and M. Potthoff, Phys. Rev. B 70 
245110 (2004). 

[27] T.A. Maier, M. Jarrell, T.C. Schulthess, P.R.C. Kent and J.B. White, Phys. Rev. Lett. 
95, 237001 (2005). 



BIBLIOGRAPHY 45 

[28] C. Huscroft, M. Jarrell, T. Maier, S. Moukouri and A.N. Tahvildarzadeh, Phys. 
Rev. Lett. 86, 139 (2001); A. Macridin, M. Jarrell, T. Maier and P.R.C. Kent, 



cond-mat/0509166j 



B. Kyung, V. Hankevych, A.-M. Dare, A.-M.S. Tremblay, Phys. Rev. Lett. 93, 147004 
(2004); B. Kyung, S.S. Kancharla, D. Senechal, A.-M. Tremblay, M. Civelli, G. Kotliar, 



|cond-mat/0502565i 



[30] C. Dahnken, M. Potthoff, E. Arrigoni and W. Hanke, cond-mat/0504618 



[31] S.S. Kancharla, M. Civelli, M. Capone, D. Senechal, G. Kotliar, A.-M.S. Tremblay, 



|cond-mat/0508205, 



[32] M. Aichhorn, E. Arrigoni, M. Potthoff and W. Hanke, |cond-mat/05 11460 



A.I. Lichtenstein and M.I. Katsnelson, cond-mat /9911320 



[34] P.W. Anderson, Science 235, 1196 (1987). 

M. Salmhofer, Commun. Math. Phys. 194, 249 (1998). 
C.J. Halboth and W. Metzner, Phys. Rev. B 61, 7364 (2000). 



[37] C. Honerkamp, M. Salmhofer, N. Furukawa and T.M. Rice, Phys. Rev. B 63, 35109 
(2001). 

[38] W.O. Putikka and M.U. Luchini, Phys. Rev. B 62, 1684 (2000). 

[39] L.P. Pryadko, S.A. Kivelson, and O. Zachar, Phys. Rev. Lett. 92, 67002 (2004). 

[40] T. Koretsune and M. Ogata, .cond-mat/0505618| 

[41] N. Bulut, D.J. Scalapino and S.R. White, Phys. Rev. B 47, R6157 (1993); ibid. Phys. 
Rev. B 50, 9623 (1994). 

[42] T.A. Maier, M. Jarrell and D.J. Scalapino, Phys. Rev. Lett. 96, 47005 (2006). 



46 BIBLIOGRAPHY 

[43] R. Blankenbecker, D.J. Scalapino and R.L. Sugar, Phys. Rev. D 24, 2278 (1981). 

[44] S.R. White, Phys. Rev. Lett. 69, 2863 (1992). 

[45] S.R. White, Phys. Rev. B 48, 10345 (1993). 

[46] J. Hirsch, Phys. Rev. Lett. 51, 1900 (1983). 

[47] E.Y. Loh, J.E. Gubernatis, R.T. Scalettar, S.R. White, D.J. Scalapino and R.L. Sugar, 
Phys. Rev. B 41, 9301 (1990). 

[48] J.E. Hirsch and R.M. Fey, Phys. Rev. Lett. 56, 2521 (1986). 

[49] M. Jarrell et al. , Phys. Rev. B 64, 195130 (2001). 

[50] D. Betts, H. Lin and J. Flynn, Can. J. Phys. 77, 353 (1999). 

[51] D.A. Huse, Phys. Rev. B 37, 2380 (1988). 

[52] J.E. Hirsch and H.Q. Lin, Phys. Rev. B 37, 5070 (1988). 

[53] In the infinite cluster limit, one expects that Pd will follow the low temperature 
Kosterlitz-Thouless behavior P^^ ~ A exp{-2B/{T - TJ"-^). 

[54] R.M. Noack, N. Bulut, D.J. Scalapino and M.G. Zacher, Phys. Rev. B 56, 7162 (1997). 

[55] O.K. Andersen, J. Phys. Chem. Solids 56, 1573 (1995); J. Low Temp. Physics 105, 285 
(1996); Phys. Rev. B 62, R16219 (2000). 

[56] E. Pavarini, I. Dasgupta, T. Saha-Dasgupta, O. Jepsen and O.K. Andersen, Phys. Rev. 
Lett. 87, 47003 (2001). 

[57] T. Maier, M. Jarrell, T. Pruschke and J. Keller, Phys. Rev. Lett. 85, 1524 (2000). 

[58] G.B. Martins, J.C. Xavier, L. Arrachea and E. Dagotto, Phys. Rev. B64, 180513 (2001). 

[59] P. Prelovsek and A. Ramsak, Phys. Rev. B 65, 174529 (2002) and cond-mat/0502044| 



BIBLIOGRAPHY 47 

[60] S.R. White and D.J. Scalapino, Phys. Rev. Lett. 91, 136403 (2003). 

[61] T. Timusk, and B. Statt, Reports of Progress in Physics 62, 61 (1999). 

[62] H. AUoul, T. Ohno and P. Mendels, Phys. Rev. Lett. 63, 1700 (1989). 

[63] Ch. Renner, B. Revaz, J.-Y. Genoud, K. Kadowaki and O. Fischer, Phys. Rev. Lett. 80, 
149 (1998). 

[64] C.C. Homes, T. Timusk, R. Liang, D.A. Bonn and W.N. Hardy, Phys. Rev. Lett. 71, 
1645 (1993). 

[65] J.W. Loram, K.A. Mirza, J.R. Cooper, W.Y. Liang and J.M. Wade, J. of Superconduc- 
tivity 7, 243 (1994). 

[66] A. DamasceUi, Z. Hussain and Z.X. Shen, Rev. Mod. Phys. 75, 474 (2004). 

[67] A. Moreo, S. Haas, A.W. Sandirk, and E. Dagotto, Phys. Rev. B 51, 12045 (1995). 

[68] D.J. Scalapino, Tr. J. Phys. 20, 560 (1996). 

[69] A. Kampf and J.R. Schrieffer, Phys. Rev. B 41,6399 (1990). 

[70] T.A. Maier, M. Jarrell and D.J. Scalapino, Phys. Rev. B 

[71] The eigenvalue for U — At on the 4-site cluster lays above the result obtained for the 
24-site cluster. As discussed in Section 3, this is expected because the 4-site cluster 

encloses just one d-wave plaquette and the embedding action does not contain pairfield 
fluctuations. 

[72] M. Chiao, R. Hill, C. Lupien, L. Taillefer, R Lambert, R. Gagon and R. Fournier, Phys. 
Rev. B62, 3554 (2000). 

[73] D.J. Van Harlinger, Rev. Mod. Phys. 67, 515 (1995). 

[74] C.C. Tsuei and J.R. Kirtley, Rev. Mod. Phys. 72, 969 (2000). 



48 BIBLIOGRAPHY 
[75] K.A. Kouznetsov et al. , Phys. Rev. Lett. 79, 3050 (1997). 

[76] H.J.H. Smilde, A. A. Golubov, Ariando, G. Rijnders, J.M. Dekkers, S. Harkema, 
D.H.A. Blank, H. Rogalla and H. Hilgenkamp, Phys. Rev. Lett. 95, 257001 (2005). 

[77] J.R. Kirtley, C.C Tsuei, Ariando, C.J.M. Verwijs, S. Harkema and H. Hilgenkamp, 

Nature Physics 

T.R Devereaux, T. Cuk, Z.-X. Shen and N. Nagaosa, Phys. Rev. Lett. 93, 117004 (2004). 
[79] CM. Varma, Phys. Rev. B 55, 14554 (1997). 

[80] P.W. Anderson, P.A. Lee, M. Randeria, T.M. Rice, N. Trivedi and F.C. Zhang, J. Phys. 
Condens. Matter 16, R755 (2004). 



P.W. Anderson, |cond-mat/0512471 



[82] D.J. Scalapino, Phys. Rep. 250, 330 (1995). 

[83] A.V. Chubukov, D. Pines and J. Schmalian, chapter in The Physics of Conventional 
and Unconventional Superconductors (edited by K.H. Bennemann and J.B. Kettersen, 
Springer- Verlag, 2002). 

[84] B. Barbiellini, O. Fischer, A.D. Kent, D.B. Mitzi and A. Kapitulnik, Physics B 194, 
1689 (1994). 

[85] J. P. Carbotte, E. Schachinger and D. Basou, Nature (London) 401, 354 (1999). 

[86] S. Pailhes, C. Ulrich, R. Fauque, V. Hinkov, Y. Sidis, A. Ivanov, CT. Lin, B. Keimer 
and P. Bourges, cond-mat /0512634 



M. Eschrig and M.R. Norman, Phys. Rev. B 67, 144503 (2003). 
T. Dahm, P.J. Hirschfeld, D.J. Scalapino and L. Zhu, Phys. Rev. B 
[89] N. Bulut and D.J. Scalapino, Phys. Rev. B 53, 5149 (1996). 



BIBLIOGRAPHY 



49 



[90] J. Brinckmann and P. A. Lee, Phys. Rev. B 65, 14502 (2002). 

[91] S. Daul, D.J. Scalapino and S.R. White, Phys. Rev. Lett. 84, 4188 (2000). 

[92] D.J. Scalapino, Physics in Canada 56, 267 (2000). 

[93] S.A. Kivelson and E. Fradkin, |3^M^at/0507459| 



[94] V.J. Emery and S.A. Kivelson, Phys. Rev. Lett. 74, 3253 (1995); E. Arrigoni, E. Fradkin 
and S.A. Kivelson, Phys. Rev. B 69, 214519 (2004). 

[95] A. Yazdani, CM. Howald, CP. Lutz, A. Kapitulnik, and D.M. Eigler, Phys. Rev. Lett. 
83, 176 (1999). 

[96] E.W. Hudson, S.H. Pan, A.K. Gupta, K.-W. Ng, and J.C Davis, Science 285, 88 (1999). 

[97] S.H. Pan, E.W. Hudson, K.M. Lang, H. Eisaki, S. Uchida, and J.C Davis, Nature 403, 
746 (2000). 

[98] J.C Bednorz and K.A. Miiller, Z. Phys. B: Condens. Matter 64, 189 (1986). 

Addendum 

This manuscript was written almost a year ago and will appear as Chapter 13 in the "Hand- 
book of High Temperature Superconductivity," edited by J.R. Schrieffer and published by 
Springer. It is narrowly focused on the results obtained from numerical studies of one par- 
ticular model, the Hubbard model. Here in this addendum [1] I would like to briefly address 
the more general questions of why the quest to understand the mechanism responsible for 
superconductivity in the high cuprates is important and where are we twenty years after 
Bednorz's and Muller's seminal discovery? As part of this, I will indicate how the numerical 
results for the Hubbard model may help to provide some answers. 

It has been suggested [2] that the answer to the first question is that we need to "do- 
mesticate the goat." That is, roughly 11,000 years ago man domesticated the wolf. During 
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the next 4000 years, different wolf species as well as various breeds of dogs were created 
and improved. But then, 7000 years ago the goat was discovered and domesticated. This 
achievement led within the next 1000 years to the domestication of sheep, cattle, chickens, 
swine, etc. and to a fundamentally different way of life. So the high Tc cuprate problem or, 
more generally, the problem of the strongly correlated electron superconductors is the goat 
and we are seeking to domesticate it. The hope is that this will lead to the development of 
new strongly correlated superconducting materials and a deeper understanding of a variety 
of existing ones. 

A short list of some of the materials [3] discussed at the M2S-HTS Dresden meeting 
which fall into the category of strongly correlated superconductors follow: 

Cuprates (hole and electron doped) 

Heavy fermions (U2 (PdPtja B, CeCoIns, PrOs4 Sbi2, PuCoGas) 

Ruthenates Sr2Ru04 

Organics kBEDT, (TMTSF)2PF6 

Cobaltates Na^CoOa (1.3 H2O) 

In addressing the second question regarding where are we with respect to understanding 
the mechanism responsible for superconductivity in the cuprates, it is useful to think back 
to a little over a decade ago. At that time a key question involved the symmetry of the gap. 
At an APS March meeting, Bertrum Batlogg was asked when we would have an experiment 
which would tell us the symmetry of the cuprate gap. His response was quick and to the 
point: We already have a number of such experiments, we just don't agree on which one 
is correct. As it turned out, we were not so far away from settling this question. At 
that time, van Harlingen and his group at the University of Illinois and Tsuei, Kirtley and 
their co-workers at IBM were carrying out phase sensitive experiments which would provide 
convincing evidence that the cuprate gap was d-wave like (for orthorhombic systems d+as) . 
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I believe that now, a decade later, the question regarding the nature of the cuprate 
pairing mechanism is at a similar stage. That is, as seen from the talks presented at this 
conference, there are clearly a number of theoretical proposals for the underlying cuprate 
pairing mechanism. Here is a short list, which I realize is incomplete (there were some 
200 theoretical talks and 700 posters covering electronic structure, many-body theories, 
phenomenology, and proposals for experiments and new materials): 

Jahn-Teller bipolarons 

Stripes (the role of inhomogeneity) 

RVB-RMFT-Gutzwiller Projected BCS 

Electron-phonon+U 

Spin-fluctuations 

Charge-fluctuations 

Electric quadropole fluctuations 

Loop current fluctuations 

dDW, dCDW 

Quantum critical point fluctuations 
Competing phases 
Pomeranchuck instabilities 
d-to-d electronic excitations 

The problem regarding the pairing mechanism is therefore reminiscent of the earlier gap 
symmetry question. It is not that there is a lack of proposals, but rather it is that we do 
not have a consensus on which one contains the appropriate description. Some of us had 
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thought that the experimental observation of a d-wave gap, which had been predicted from 
the analysis of Hubbard and t-J models, would have provided "convincing" evidence that 
the pairing arose from a spin mediated interaction. But as we now know, a dx^_y2 gap is not 
a sufficiently unique signature. There are a number of possible pairing mechanisms such as 
Big phonons, fluctuating current loops, as well as others which may drive d-wave pairing. 
In addition, there is the suggestion that the pairing interaction is amplified by intrinsic 
inhomogeneities or stripes. There is also the question of the possible role of amplification 
that may be obtained near a quantum critical point (QCP). Although, I believe that this 
latter QCP scenario can be viewed as tuning the parameters to optimize a given pairing 
mechanism, much as what was done when Tc was increased in the electron-phonon systems 
by varying the composition in order to be near a lattice instability. In any event, the question 
remains, how will we know? Let's consider some possibilities. 

Perhaps it will be shown that there is another ordered phase associated with the pseudo- 
gap regime such as a d-DW or a d-CDW, a time reversal breaking current-loop phase (recent 
neutron scattering and Kerr effect experiments) , a non-superconducting xy-like phase (non- 
analytic dependence of the magnetization as H ^ O at temperatures above Tc), an exotic 
spin-liquid phase, . . . There were in fact both theoretical and experimental talks on such 
phases, and definitive evidence for a new ordered phase would certainly narrow the theo- 
retical possibilities. Here however it is interesting to note that in spite of the evidence for 
a striped phase in Lai.35Sro.25Ndo.4Cu04 and increasing evidence for stripe-like fiuctuations 
in the cuprates, the question of whether stripes, or other intrinsic inhomogeneities, play an 
essential role in the high Tc pairing mechanism remains open. We did hear about ARPES 
experiments on La2-j;Srj;Cu04 which found that the gap peaks near x — 1/8 where Tc dips. 
Further experiments on dynamic inhomogeneity-induced pairing and the role of mesoscale 
structures in the high Tc problem are needed. 

Perhaps we will know when we have a further understanding of the clues coming from 
the electronic structure calculations regarding the variation of Tc with the chemical structure 
and the effective band parameters such as the next near neighbor hopping t'. Perhaps the 
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insight we need will come from an understanding of the local density of states measured by 
STM and its relation to nearby local structural and chemical changes. Perhaps the ARPES, 
INS and STM studies will show that phonons play an important role. If this is the case, 
it seems likely that it will not be in the traditional way they do in the low Tc metals, but 
rather in combination with the stripes or mesoscale structures. 

It may be, as suggested at this meeting, that one can identify a quantum critical point 
(QCP) which is associated with superconductivity and then from the nature of the QCP, 

identify the fluctuations which mediate the pairing. Of course if the interplay between the 
various phases is sufficiently strong, one may have to reconsider the meaning of mechanism. 

One route that has been discussed is that of identifying the cuprate pairing mechanism by 
showing the similarity of the cuprates to other classes of materials which exhibit unconven- 
tional superconductivity such as the actinide metals (PuCoGas, PuRhGas) and the heavy 
fermion systems (GeRIns, R=Go,Rh,Ir). Perhaps there is a common mechanism which is 
tunned as one moves from one system to another. The candidate mechanism discussed at 
this meeting was spin-fluctuations. 

Another approach is based on using numerical techniques (or possibly experimental cold 
atom analogues) to study particular models such as the Hubbard model. Here the idea is 
to determine first whether a given model exhibits the phenomena seen in the cuprates and 
then, if it does, determine the nature of the pairing interaction in the model. (This is the 
approach taken in this chapter.) 

Perhaps the best approach will be to follow the path taken for the traditional low temper- 
ature superconductors where the structure of the pairing interaction was determined from 
the tunneling dl/dV characteristic. Here of course one knew how to get A (a;) from dl/dV 
and had the Ehashberg theory to determine the interaction from A (a;). For the cuprates 
one will likely need a combination of numerical and analytic approaches along with ARPES, 
INS, STM and conductivity data to carry out such a program. 

For the cuprates, the k dependence of the gap near the fermi surface is known from 
a variety of experiments. A detailed map for YBGO is now available from 7r-junction in- 



54 



BIBLIOGRAPHY 



terference measurements which find an anisotropic d+as wave form. If one had a simple 
(cos kx — cos ky) form over the entire Brillouin zone, one would know that the pairing arose 
from an attractive near-neighbor interaction. While it is likely somewhat more extended and, 
for the case of an orthorhombic crystal, anisotropic, the simple d-wave form is a reasonable 
starting point. What is needed are experimental determinations of the frequency dependence 
of the gap. Actually one would like to determine both the renormalization parameter Z{k, lu) 
and A{k,u!) — (f){k,uj)/Z{k,uj) with (f){k,uj) the gap parameter. Studies aimed at extracting 
the frequency dependence of the gap and the underlying interaction using conductivity (t{u}) 
as well as INS and ARPES data were reported. One can expect further progress in this 
direction as such data becomes available on the same crystals. 

So at the present time, the question of the mechanism responsible for pairing in the high 
Tc cuprates remains open. However, it is clear from the work presented at this conference 
that we are moving closer to an understanding. 
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